Mobile Footprinting

1. Setup Environment

require(ggplot2)
require(summarytools)
require(cowplot)
require(caret)
require(corrplot)
require(RColorBrewer)
require(vembedr)
require(Rmisc)
require(varian)

package ‘rstan’ was built under R version 4.0.2Failed with error:  ‘package ‘rstan’ could not be loaded’

require(patchwork)
require(plotly)
require(Metrics)
require(dplyr)
require(ggpubr)
require(mosaic)
require(openxlsx)
require(visreg)
require(factoextra)
require(rstatix)
require(gridExtra)
require(colorspace)
require(grid)
require(data.table)
require(psych)
source('~/Documents/GitHub/SmartphoneSensorPipeline/Extra/plotting_functions.R')

project_path = "~/Documents/xia_gps/"
data_path = file.path(project_path,"beiwe_output_043020")
gps_df_path = file.path(data_path,"Processed_Data/Group/feature_matrix.txt")

2.Prepare GPS data

gps_df = read.table(gps_df_path,header = T, dec = ",")[,c(1,2,97:111)]
gps_df[,3:17] = apply(gps_df[,3:17],2,as.numeric)

# loop through each subj to remove 1st and last days of gps data
    gps_df_clean = data.frame() #initiate a df
    for (subj in unique(gps_df$IID)){ #loop through each subj
      gps_df_subj <- subset(gps_df, IID == subj) #get gps_df per subject
      gps_df_subj <- gps_df_subj[2:(dim(gps_df_subj)[1]-1),] #remove the 1st and last days
      gps_df_clean <- rbind(gps_df_clean,gps_df_subj) #combine all subjs
    }

    # exclude days with too much excessiveness
    sensitivity_cutoff = 1440 # this controls the cutoff threshold
    gps_df_clean2 = subset(gps_df_clean, MinsMissing < sensitivity_cutoff)

set.seed(510)
make_subj_seq = function(gps_df_clean2,part_times) {
  subj_seq = list()
  for (subj in unique(gps_df_clean2$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    if (dim(subj_data)[1] >5) {
    subj_seq[[subj]] <-createDataPartition(subj_data$IID,times = part_times, p =0.5)
    }
  }
  return(subj_seq)
}
part_times = 1000
subj_seq = make_subj_seq(gps_df_clean2, part_times)

# an example of subj 1, and first half
example_data = gps_df_clean2[subj_seq$`16xv6ko1`$Resample0001,3:17]
gps_cor = rquery.cormat(example_data, type = "full")

gps_clean2_feature = make_feature_matrix(gps_df_clean2, subj_seq, 3:17 )
subject_seq = names(gps_clean2_feature$subj_mat_1)

gps_wide_matrix = data.frame()
subj_days = gps_df_clean2 %>% group_by(IID) %>% dplyr::tally(name = "Days Collected")
for (subj in arrange(subj_days,`Days Collected`)$IID){
  #if ((subj %in% subj_days$IID[which(subj_days$`Days Collected`<=10)]) == T) {
    for (part in 1:5){
      half_1 = gps_clean2_feature$subj_mat_1[[subj]][[part]]$cor
      half_2 = gps_clean2_feature$subj_mat_2[[subj]][[part]]$cor
      half_1_half_2 = c(half_1,half_2)
      gps_wide_matrix = rbind(gps_wide_matrix,half_1)
      gps_wide_matrix = rbind(gps_wide_matrix,half_2)
      }
    #}
}

gps_wide_matrix = t(gps_wide_matrix)
gps_corplot = rquery.cormat(gps_wide_matrix, type = "full",graph=FALSE)
gps_corplot$subj = arrange(subj_days,`Days Collected`)$IID

levelplot(gps_corplot$r,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Pearson correlation", main=list(label='GPS Feature Similarity'),xlab="",ylab="")

3. Prepare accelerometer data

subj_list = unique(gps_df$IID)
subj_acc_data_dim = dim(rbindlist(lapply(subj_list, function(subj) {subj_data =  readRDS(file.path(data_path,"Results/Group/accelerometer",subj,"accelerometer_ft.rds")); subj_data[,c(2:8,11)]})))

subj_acc_data_clean = as.data.frame(rbindlist(lapply(subj_list, function(subj) { subj_data =  readRDS(file.path(data_path,"Results/Group/accelerometer",subj,"accelerometer_ft.rds")); subj_data[2:(dim(subj_data)[1]-1),c(2:8,11)]}))) #also removes 1st and last day

example_acc_data = subj_acc_data_clean[subj_seq$`16xv6ko1`$Resample0001,1:7]

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls

acc_cor = rquery.cormat(example_acc_data, type = "full")

subj_acc_data_clean$IID = subj_acc_data_clean$subject
set.seed(510)
acc_subj_seq = make_subj_seq(subj_acc_data_clean, part_times)
acc_clean_feature = make_feature_matrix(subj_acc_data_clean, acc_subj_seq,  c(1:7) )

acc_wide_matrix = data.frame()
subj_acc_days = subj_acc_data_clean %>% group_by(IID) %>% dplyr::tally(name = "acc_days")
for (subj in arrange(subj_acc_days,`acc_days`)$IID){
    for (part in 1:5){
      half_1 = acc_clean_feature$subj_mat_1[[subj]][[part]]$cor
      half_2 = acc_clean_feature$subj_mat_2[[subj]][[part]]$cor
      half_1_half_2 = c(half_1,half_2)
      acc_wide_matrix = rbind(acc_wide_matrix,half_1)
      acc_wide_matrix = rbind(acc_wide_matrix,half_2)
      }
}

acc_wide_matrix = t(acc_wide_matrix)
acc_corplot = rquery.cormat(acc_wide_matrix, type = "full",graph=FALSE)
acc_corplot$subj = arrange(subj_acc_days,`acc_days`)$IID

levelplot(acc_corplot$r,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Pearson correlation", main=list(label='Acc Feature Similarity'),xlab="",ylab="")

4. Merge GPS and accelerometer data

gps_acc_combined_feat = list()
for (subj in subj_list){
  for (time in 1:part_times){
    gps_feat1 = gps_clean2_feature$subj_mat_1[[subj]][[time]]$cor
    acc_feat1 = acc_clean_feature$subj_mat_1[[subj]][[time]]$cor
    
    gps_feat2 = gps_clean2_feature$subj_mat_2[[subj]][[time]]$cor
    acc_feat2 = acc_clean_feature$subj_mat_2[[subj]][[time]]$cor

    gps_acc_combined_feat$subj_mat_1[[subj]][[time]] = c(gps_feat1, acc_feat1)
    gps_acc_combined_feat$subj_mat_2[[subj]][[time]] = c(gps_feat2, acc_feat2)
  }
}

5. GPS and accelerometer based Individual Identification

match_combined_gps_accfeat = calc_match_vector(gps_acc_combined_feat$subj_mat_1, gps_acc_combined_feat$subj_mat_2, "cor")
acc_time_cb_accgps  = calc_acc_time(match_combined_gps_accfeat, "max")
acc_subj_cb_accgps  = calc_acc_subj(match_combined_gps_accfeat, "max")

p = hist_chx(acc_time_cb_accgps, bins = 14, title = paste("Individual Identificaiton Accuracy \n (GPS+ Accel)"), xaxis = "Prediction Accuracy", yaxis = "Count")
ggplotly(p)

perm_time = 10
gps_df_perm = gps_df_clean2
acc_df_perm = subj_acc_data_clean

perm_acc_gps_time = list()
perm_acc_gps_subj = list()

for (i in 1:perm_time) {
  perm_part_times = 1
  print(paste("processing ...", i,"..."))
  gps_df_perm$IID = sample(gps_df_perm$IID)
  acc_df_perm$IID = sample(acc_df_perm$IID)
  
  perm_subj_seq = make_subj_seq(gps_df_perm,part_times = perm_part_times)
  perm_acc_subj_seq = make_subj_seq(acc_df_perm,part_times = perm_part_times)
  
  perm_gps_cor = make_feature_matrix(gps_df_perm,perm_subj_seq,3:17)
  perm_acc_cor = make_feature_matrix(acc_df_perm,perm_acc_subj_seq,1:7)
  
  subj_mat_1 = list()
  subj_mat_2 = list()
  
    for (subj in subj_list){
        subj_mat_1[[subj]]$Resample1 = as.numeric(c(perm_gps_cor$subj_mat_1[[subj]][[1]]$cor, perm_acc_cor$subj_mat_1[[subj]][[1]]$cor))
        subj_mat_2[[subj]]$Resample1 = as.numeric(c(perm_gps_cor$subj_mat_2[[subj]][[1]]$cor, perm_acc_cor$subj_mat_2[[subj]][[1]]$cor))
    }
  perm_gps_acc_mats = list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2)
  
  perm_match = calc_match_vector(perm_gps_acc_mats$subj_mat_1,perm_gps_acc_mats$subj_mat_2, "cor")
  perm_acc_gps_time[[i]] = calc_acc_time(perm_match, "max")
  perm_acc_gps_subj[[i]] = calc_acc_subj(perm_match, "max")
}


perm_acc_gps_time = unlist(perm_acc_gps_time)
perm_acc_gps_subj = unlist(perm_acc_gps_subj)

df = data.frame(val = perm_acc_gps_time)

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls

p = ggplot(df, aes(x=val)) + 
    geom_histogram(bins = 7, fill = "#9ECAE1", size = 2) + 
    theme_cowplot() +
    labs(title=paste("GPS+Accel: \n",part_times,"Permutations"),
         x  = "Individual Identificaiton Accuracy", y = "Count") + 
  #geom_density(alpha=.2) +
    theme(plot.title = element_text(hjust = 0.5))
ggplotly(p)

perm_subj = sapply(subject_seq, function(subj) sum(perm_acc_gps_subj[which(names(perm_acc_gps_subj) == subj)])/1000)

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls

subj_df <- data.frame(x=names(subj_seq))
subj_df$y = acc_subj_cb_accgps
subj_df$y_perm = perm_subj
subj_df = subj_df[order(subj_df$y),]
p_subj_cor_perm = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
  geom_point(aes(y = y), color = "#F69274", size = 2) + 
  geom_point(aes(y = y_perm), color = "#9ECAE1", size =2 ) +
  theme_cowplot() + 
  labs(title = paste("GPS+Accel \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
       x = "Subjects", y = "Individual Identification Accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(p_subj_cor_perm)

6. GPS only based Individual Identification

gps_match_cor = calc_match_cor(gps_clean2_feature$subj_mat_1,gps_clean2_feature$subj_mat_2)
gps_acc_time = calc_acc_time(gps_match_cor, "max")
gps_acc_subj = calc_acc_subj(gps_match_cor)

gps_df_perm = gps_df_clean2
perm_time = 1000
perm_acc_time = list()
perm_acc_subj = list()
for (i in 1:perm_time) {
  perm_part_times = 1
  print(paste("processing ...", i,"..."))
  gps_df_perm$IID = sample(gps_df_perm$IID)
  perm_subj_seq = make_subj_seq(gps_df_perm, part_times = perm_part_times)
  perm_gps = make_feature_matrix(gps_df_perm,perm_subj_seq,3:17)
  perm_mat_1 = perm_gps$subj_mat_1
  perm_mat_2 = perm_gps$subj_mat_2
  perm_match_cor = calc_match_cor(perm_mat_1,perm_mat_2)
  perm_acc_time[[i]] = calc_acc_time(perm_match_cor)
  perm_acc_subj[[i]] = calc_acc_subj(perm_match_cor)
}

p_time_cor = hist_chx(acc_time, bins = 17, title = paste("GPS: \n Individual identification accuracy"), xaxis = "Individual identification accuracy", yaxis = "Count")
ggplotly(p_time_cor)

#perm_acc_time_all = unlist(perm_acc_time)

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls

q_time_cor = hist_chx(perm_acc_time_all, bins = 7, title = paste("GPS accuracy across",length(perm_acc_time_all),"Permutations"), xaxis = "Individual Identificaiton Accuracy", yaxis = "Count")
ggplotly(q_time_cor)

perm_acc_subj = unlist(perm_acc_subj)
perm_subj = sapply(subject_seq, function(subj) sum(perm_acc_subj[which(names(perm_acc_subj) == subj)])/1000)

subj_df <- data.frame(x=names(subj_seq))
subj_df$y = gps_acc_subj
subj_df$y_perm = perm_subj
subj_df = subj_df[order(subj_df$y),]
p_subj_cor_perm = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
  geom_point(aes(y = y, col = "subject data")) + 
  geom_point(aes(y = y_perm, col = "permutation")) +
  theme_cowplot() + 
  labs(title = paste("GPS \n Subject Level Accuracy Across",perm_time,"permutations"), 
       x = "Subjects", y = "Individual identification accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(p_subj_cor_perm)

7. Accelerometer only based Individual Identification

acc_match_cor = calc_match_cor(acc_clean_feature$subj_mat_1,acc_clean_feature$subj_mat_2)
acc_acc_time = calc_acc_time(acc_match_cor, "max")
acc_acc_subj = calc_acc_subj(acc_match_cor)

acc_df_perm = subj_acc_data_clean
perm_time = 1000
perm_acc_acc_time = list()
perm_acc_acc_subj = list()
for (i in 1:perm_time) {
  perm_part_times = 1
  print(paste("processing ...", i,"..."))
  acc_df_perm$IID = sample(acc_df_perm$IID)
  perm_subj_seq = make_subj_seq(acc_df_perm, part_times = perm_part_times)
  perm_acc = make_feature_matrix(acc_df_perm,perm_subj_seq,1:7)
  perm_mat_1 = perm_acc$subj_mat_1
  perm_mat_2 = perm_acc$subj_mat_2
  perm_match_cor = calc_match_cor(perm_mat_1,perm_mat_2)
  perm_acc_acc_time[[i]] = calc_acc_time(perm_match_cor)
  perm_acc_acc_subj[[i]] = calc_acc_subj(perm_match_cor)
}

df = data.frame(val = acc_acc_time)
p = ggplot(df, aes(x=val)) + 
    geom_histogram(bins = 15, fill = "#F69274", size = 2) + 
    theme_cowplot() +
    labs(title=paste("Accelerometer Features: \n Average Accuracy Across",part_times,"Data Partitions"),
         x  = "Mobility Footprint Individual Identification Accuracy", y = "Count") + 
  #geom_density(alpha=.2) +
    theme(plot.title = element_text(hjust = 0.5))
ggplotly(p)

df = data.frame(val = perm_acc_acc_time)
p = ggplot(df, aes(x=val)) + 
    geom_histogram(bins = 5, fill = "#9ECAE1", size = 2) + 
    theme_cowplot() +
    labs(title=paste("Accelerometer Features: \n Average Accuracy Across",perm_time,"Permutations"),
         x  = "Mobility Footprint Individual Identification Accuracy", y = "Count") + 
  #geom_density(alpha=.2) +
    theme(plot.title = element_text(hjust = 0.5))
ggplotly(p)

perm_acc_acc_subj = unlist(perm_acc_acc_subj)

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls

perm_subj = sapply(subject_seq, function(subj) sum(perm_acc_acc_subj[which(names(perm_acc_acc_subj) == subj)])/1000)

subj_df <- data.frame(x=names(subj_seq))
subj_df$y = acc_acc_subj
subj_df$y_perm = perm_subj
subj_df = subj_df[order(subj_df$y),]
p_subj_cor_perm = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
  geom_point(aes(y = y, col = "subject data")) + 
  geom_point(aes(y = y_perm, col = "permutation")) +
  theme_cowplot() + 
  labs(title = paste("Accelerometer \n Subject Level Accuracy Across",perm_time,"permutations"), 
       x = "Subjects", y = "Individual identification accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(p_subj_cor_perm)

df_accgps_time = data.frame(GPS = gps_acc_time, Accelerometer = acc_acc_time, `GPS and Acceleromter` = acc_time_cb_accgps)
df_accgps_time_long = gather(df_accgps_time, feature, accuracy, 1:3, factor_key=TRUE)
mu = ddply(df_accgps_time_long, "feature", summarise, grp.mean=mean(accuracy))
p_accgps_feature = ggplot(df_accgps_time_long, aes(x=accuracy,color = feature, fill = feature)) + 
    geom_density(size = 1, alpha = 0) + 
    #geom_histogram(aes(y=..density..), alpha=0.2, 
    #            position="identity") + 
    geom_vline(data=mu, aes(xintercept=grp.mean, color=feature),
             linetype="dashed") +
    labs(y = "Density", x = "Individual Identitication Accuracy") + 
    theme_cowplot() +
    theme(legend.position="top")
p_accgps_feature

8. Data quantity vs. quality associations

gps_days = gps_df_clean2 %>% group_by(IID) %>% dplyr::tally(name = "gps_days")

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls

acc_days = subj_acc_data_clean %>% group_by(IID) %>% dplyr::tally(name = "acc_days")
data_quant_plots = list()
data_quant_plots$gps = conf_scatter_plot(gps_days, gps_acc_subj, "gps_days")
data_quant_plots$acc = conf_scatter_plot(acc_days, acc_acc_subj, "acc_days")
Reduce(`+`,data_quant_plots)

gps_quality = gps_df_clean2 %>% group_by(IID) %>% dplyr::summarize(MinsMissing = mean(MinsMissing,na.rm=T))
#acc_quality = sapply(subj_list,get_data_mean)
#acc_quality$IID = names(acc_data_quality)
#acc_quality$datapoints = as.numeric(value(unlist(acc_data_quality)[1:41]))

data_qual_plots = list()
data_qual_plots$gps = conf_scatter_plot(gps_quality, gps_acc_subj, "MinsMissing")
data_qual_plots$acc = conf_scatter_plot(acc_quality, acc_acc_subj, "datapoints")

Reduce(`+`,data_qual_plots)

9. Developmental effects and sex differences

ids = read.xlsx(file.path(project_path,"data/clinical_data/subjecttracker_4.xlsx"))[1:41,c(1,3)]
demo_tem = read.csv(file.path(project_path,"data/clinical_data/self_report_itemwise.csv"), colClasses = c("NULL",NA,"NULL","NULL",NA,NA,rep("NULL",363)))

demo_beiwe = inner_join(ids,demo_tem, by = c("BBLID" = "bblid"))
acc_subj_df = data.frame(beiweID = names(acc_subj), gps = gps_acc_subj, acc = acc_acc_subj, accgps = acc_subj_cb_accgps)
acc_demo_df = inner_join(acc_subj_df,demo_beiwe, by = "beiweID")
acc_demo_df = inner_join(acc_demo_df,gps_days, by = c("beiweID" = "IID"))
acc_demo_df = inner_join(acc_demo_df,acc_days, by = c("beiweID" = "IID"))
acc_demo_df$admin_sex = as.factor(acc_demo_df$admin_sex)

agesex_gps_fit = gam(gps ~ s(admin_age) + admin_sex + gps_days, data = acc_demo_df, method = "REML")

p_gps = list()
p_gps$sex = visreg(agesex_gps_fit, "admin_sex", gg = T, line=list(col="#3576b5",size = 0.5, alpha = 0.25), 
       points=list(size=2, pch=19), fill= list(fill=c("#b3d3f2"), alpha = 0.25)) + 
  xlab("Sex") + ylab("Individual Footprint Distinctiveness \n (GPS)") +  theme_cowplot()

p_gps$age = visreg(agesex_gps_fit, "admin_age", gg = T, line=list(col="#3576b5",size = 0.5, alpha = 0.25), 
       points=list(size=2, pch=19), fill= list(fill=c("#b3d3f2"), alpha = 1)) + 
  xlab("Age") + ylab("Individual Footprint Distinctiveness \n (GPS)") +  theme_cowplot()

Reduce(`+`,p_gps)

agesex_acc_fit = gam(acc ~ s(admin_age) + admin_sex + acc_days , data = acc_demo_df, method = "REML")

p_acc = list()

p_acc$age = visreg(agesex_acc_fit, "admin_age", gg = T, line=list(col="#3576b5",size = 0.5, alpha = 0.25), 
       points=list(size=2, pch=19), fill= list(fill=c("#b3d3f2"), alpha = 0.25)) + 
  xlab("Age") + ylab("Individual Footprint Distinctiveness \n (Accel)") +  theme_cowplot()

p_acc$sex = visreg(agesex_acc_fit, "admin_sex", gg = T, line=list(col="#3576b5",size = 0.5, alpha = 0.25), 
       points=list(size=2, pch=19), fill= list(fill=c("#b3d3f2"), alpha = 1)) + 
  xlab("Sex") + ylab("Individual Footprint Distinctiveness \n (Accel)") +  theme_cowplot()

Reduce(`+`,p_acc)

10. Mood and sleep variability

11. Functional connectivity